Intro

Fit density (also referred to as CPUE) model with environmental predictors and use that to calculate weighted mean dissolved oxygen, temperature and depth of Baltic cod

# Load libraries, install if needed
library(tidyverse); theme_set(theme_classic(base_size = 10))
library(readxl)
library(tidylog)
library(RCurl)
library(viridis)
library(RColorBrewer)
library(patchwork)
library(janitor)
library(icesDatras)
library(mapdata)
library(patchwork)
library(rgdal)
library(raster)
library(sf)
library(rgeos)
library(chron)
library(lattice)
library(ncdf4)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(EnvStats)
library(qwraps2)
library(bayesplot)
library(tidybayes)
library(brms)
#remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB) # To load entire cache in interactive r session, do: qwraps2::lazyload_cache_dir(path = "R/analysis/cod_flounder_diets_spatial_cache/html")

For maps

# Specify map ranges
ymin = 55; ymax = 58; xmin = 14; xmax = 20

map_data <- rnaturalearth::ne_countries(
  scale = "medium",
  returnclass = "sf", continent = "europe")

# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
swe_coast <- suppressWarnings(suppressMessages(
  st_crop(map_data,
          c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))

# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)

ggplot(swe_coast_proj) + geom_sf() 


# Define plotting theme for main plot
theme_plot <- function(base_size = 10, base_family = "") {
  theme_light(base_size = 10, base_family = "") +
    theme(
      axis.text.x = element_text(angle = 90),
      axis.text = element_text(size = 8),
      legend.text = element_text(size = 8),
      legend.title = element_text(size = 8),
      legend.position = "bottom",
      legend.key.height = unit(0.2, "cm"),
      legend.margin = margin(0, 0, 0, 0),
      legend.box.margin = margin(-5, -5, -5, -5),
      strip.text = element_text(size = 8, colour = 'black', margin = margin()),
      strip.background = element_rect(fill = "grey90")
      )
}

# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 10, base_family = "") {
  theme_light(base_size = 10, base_family = "") +
    theme(
        axis.text.x = element_text(angle = 90),
        axis.text = element_text(size = 6),
        strip.text = element_text(size = 8, colour = 'black', margin = margin()),
        strip.background = element_rect(fill = "grey90")
      )
}

# Make default base map plot
plot_map_raster <- 
ggplot(swe_coast_proj) + 
  geom_sf(size = 0.3) +
  labs(x = "Longitude", y = "Latitude") +
  theme_facet_map(base_size = 14)

# Function to convert from lat lon to UTM
LongLatToUTM <- function(x, y, zone){
  xy <- data.frame(ID = 1:length(x), X = x, Y = y)
  coordinates(xy) <- c("X", "Y")
  proj4string(xy) <- CRS("+proj=longlat +datum=WGS84")  ## for example
  res <- spTransform(xy, CRS(paste("+proj=utm +zone=",zone," ellps=WGS84",sep='')))
  return(as.data.frame(res))
}

theme_clean <- function() {
  theme_minimal(base_family = "Barlow Semi Condensed") +
    theme(panel.grid.minor = element_blank(),
          plot.title = element_text(family = "BarlowSemiCondensed-Bold"),
          axis.title = element_text(family = "BarlowSemiCondensed-Medium"),
          strip.text = element_text(family = "BarlowSemiCondensed-Bold",
                                    size = rel(1), hjust = 0),
          strip.background = element_rect(fill = "grey80", color = NA))
}

Read and plot data

# These data are for total and prey specific stomach models
cod <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/cod_diet_analysis.csv") %>% dplyr::select(-X1)
#> Warning: Missing column names filled in: 'X1' [1]
#> 
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#>   .default = col_double(),
#>   unique_pred_id = col_character(),
#>   time_period = col_character(),
#>   quarter_fact = col_character(),
#>   predator_spec = col_character(),
#>   ices_rect = col_character(),
#>   cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.

cod <- cod %>%
  mutate(year = as.integer(year),
         quarter = as.factor(quarter),
         depth2_sc = depth - mean(depth),
         saduria_entomon_per_mass = saduria_entomon_tot/pred_weight_g,
         tot_prey_biom_per_mass = tot_prey_biom/pred_weight_g,
         depth_sc = (depth - mean(depth)) / sd(depth)) %>% 
  filter(year > 2014) %>%
  filter(!quarter == 2) %>% 
  drop_na(predfle_density_sc, predcod_density_sc) %>% 
  droplevels()
#> mutate: converted 'year' from double to integer (0 new NA)
#>         converted 'quarter' from double to factor (0 new NA)
#>         new variable 'depth2_sc' (double) with 106 unique values and 0% NA
#>         new variable 'saduria_entomon_per_mass' (double) with 1,709 unique values and 0% NA
#>         new variable 'tot_prey_biom_per_mass' (double) with 7,202 unique values and 0% NA
#>         new variable 'depth_sc' (double) with 106 unique values and 0% NA
#> filter: removed 4,968 rows (57%), 3,703 rows remaining
#> filter: removed 91 rows (2%), 3,612 rows remaining
#> drop_na: no rows removed

fle <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/fle_diet_analysis.csv") %>% dplyr::select(-X1)
#> Warning: Missing column names filled in: 'X1' [1]
#> 
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#>   .default = col_double(),
#>   unique_pred_id = col_character(),
#>   predator_spec = col_character(),
#>   ices_rect = col_character(),
#>   cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.

fle <- fle %>%
  mutate(year = as.integer(year),
         quarter = as.factor(quarter),
         depth2_sc = depth - mean(depth),
         saduria_entomon_per_mass = saduria_entomon_tot/pred_weight_g,
         tot_prey_biom_per_mass = tot_prey_biom/pred_weight_g,
         depth_sc = (depth - mean(depth)) / sd(depth)) %>% 
  filter(!quarter == 2) %>% 
  drop_na(predfle_density_sc, predcod_density_sc) %>% 
  droplevels()
#> mutate: converted 'year' from double to integer (0 new NA)
#>         converted 'quarter' from double to factor (0 new NA)
#>         new variable 'depth2_sc' (double) with 81 unique values and 0% NA
#>         new variable 'saduria_entomon_per_mass' (double) with 608 unique values and 0% NA
#>         new variable 'tot_prey_biom_per_mass' (double) with 1,735 unique values and 0% NA
#>         new variable 'depth_sc' (double) with 81 unique values and 0% NA
#> filter: removed 58 rows (3%), 2,190 rows remaining
#> drop_na: no rows removed

# These data are for the diet analysis (all prey groups)
cod_prey <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/pca_cod_diet_analysis.csv") %>% dplyr::select(-X1) %>% mutate(species = "COD") %>% filter(year > 2014) %>% filter(!quarter == 2) %>% droplevels()
#> Warning: Missing column names filled in: 'X1' [1]
#> 
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#>   .default = col_double(),
#>   unique_pred_id = col_character(),
#>   predator_spec = col_character(),
#>   ices_rect = col_character(),
#>   cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> filter: removed 4,968 rows (57%), 3,703 rows remaining
#> filter: removed 91 rows (2%), 3,612 rows remaining

fle_prey <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/pca_fle_diet_analysis.csv") %>% dplyr::select(-X1) %>% mutate(species = "FLE") %>% filter(!quarter == 2) %>% droplevels()
#> Warning: Missing column names filled in: 'X1' [1]
#> 
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#>   .default = col_double(),
#>   unique_pred_id = col_character(),
#>   predator_spec = col_character(),
#>   ices_rect = col_character(),
#>   cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> filter: removed 58 rows (3%), 2,190 rows remaining

Read the prediction grids

# And now read in pred_grid2 which has oxygen values at location and time and depth:
pred_grid2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid2.csv") ### CHANGE TO OUR NEW PRED GRID

# Clean
pred_grid2 <- pred_grid2 %>%
  mutate(year = as.integer(year)) %>% 
  drop_na(depth)

# Add ices_rect
pred_grid2$ices_rect <- ices.rect2(pred_grid2$lon, pred_grid2$lat) 

# pred_grid2_q1 <- pred_grid2 %>% mutate(quarter = factor(1))
# pred_grid2_q4 <- pred_grid2 %>% mutate(quarter = factor(4))

Calculate Schoener’s overlap index for diet

# Plot data in space
plot_map_raster +
  geom_point(data = fle_prey, aes(x = X * 1000, y = Y * 1000, color = "FLE"), alpha = 0.5) +
  geom_point(data = cod_prey, aes(x = X * 1000, y = Y * 1000, color = "COD"), alpha = 0.5) +
  theme_plot()


cod_prey <- cod_prey %>% filter(lat < 58)
#> filter: removed 19 rows (1%), 3,593 rows remaining
fle_prey <- fle_prey %>% filter(lat < 58)
#> filter: removed 10 rows (<1%), 2,180 rows remaining

# Reformat data to calculate Schoeners overlap index
# colnames(fle_prey)
fle_prey_long <- fle_prey %>%
  pivot_longer(14:29) %>% 
  group_by(name, year, quarter, ices_rect) %>% # ices rect?
  summarise(fle_diet = sum(value)) %>% 
  arrange(name, year, ices_rect) %>%
  ungroup() %>% 
  group_by(year, quarter, ices_rect) %>%
  mutate(fle_diet_tot = sum(fle_diet),
         fle_diet_prop = fle_diet / fle_diet_tot) %>% 
  ungroup()
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 2180x30, now 34880x16]
#> group_by: 4 grouping variables (name, year, quarter, ices_rect)
#> summarise: now 1,296 rows and 5 columns, 3 group variables remaining (name, year, quarter)
#> ungroup: no grouping variables
#> group_by: 3 grouping variables (year, quarter, ices_rect)
#> mutate (grouped): new variable 'fle_diet_tot' (double) with 81 unique values and 0% NA
#>                   new variable 'fle_diet_prop' (double) with 445 unique values and 1% NA
#> ungroup: no grouping variables

fle_prey_wide <- fle_prey_long %>%
  pivot_wider(names_from = name, values_from = fle_diet_prop, values_fill = 0)
#> pivot_wider: reorganized (name, fle_diet_prop) into (amphipoda_tot, bivalvia_tot, clupea_harengus_tot, clupeidae_tot, gadiformes_tot, …) [was 1296x7, now 526x21]

cod_prey_long <- cod_prey %>%
  pivot_longer(14:29) %>% 
  group_by(name, year, quarter, ices_rect) %>%
  summarise(cod_diet = sum(value)) %>% 
  arrange(name, year, ices_rect) %>%
  ungroup() %>% 
  group_by(year, quarter, ices_rect) %>%
  mutate(cod_diet_tot = sum(cod_diet),
         cod_diet_prop = cod_diet / cod_diet_tot) %>% 
  ungroup()
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 3593x30, now 57488x16]
#> group_by: 4 grouping variables (name, year, quarter, ices_rect)
#> summarise: now 1,696 rows and 5 columns, 3 group variables remaining (name, year, quarter)
#> ungroup: no grouping variables
#> group_by: 3 grouping variables (year, quarter, ices_rect)
#> mutate (grouped): new variable 'cod_diet_tot' (double) with 106 unique values and 0% NA
#>                   new variable 'cod_diet_prop' (double) with 784 unique values and 0% NA
#> ungroup: no grouping variables

cod_prey_wide <- cod_prey_long %>%
  pivot_wider(names_from = name, values_from = cod_diet_prop, values_fill = 0)
#> pivot_wider: reorganized (name, cod_diet_prop) into (amphipoda_tot, bivalvia_tot, clupea_harengus_tot, clupeidae_tot, gadiformes_tot, …) [was 1696x7, now 892x21]

# Calculate Schoener index
schoener <- left_join(cod_prey_long, fle_prey_long) %>% 
  drop_na(name) %>% 
  drop_na(fle_diet_prop) %>% 
  drop_na(cod_diet_prop) %>% 
  group_by(year, quarter, ices_rect) %>%
  summarise(schoener = 1 - 0.5*(sum(abs(fle_diet_prop - cod_diet_prop)))) %>% 
  ungroup()
#> Joining, by = c("name", "year", "quarter", "ices_rect")
#> left_join: added 3 columns (fle_diet, fle_diet_tot, fle_diet_prop)
#>            > rows only in x     432
#>            > rows only in y  (   32)
#>            > matched rows     1,264
#>            >                 =======
#>            > rows total       1,696
#> drop_na: no rows removed
#> drop_na: removed 448 rows (26%), 1,248 rows remaining
#> drop_na: no rows removed
#> group_by: 3 grouping variables (year, quarter, ices_rect)
#> summarise: now 78 rows and 4 columns, 2 group variables remaining (year, quarter)
#> ungroup: no grouping variables

# Calculate Levin niche index
levin <- left_join(cod_prey_long, fle_prey_long) %>% 
  drop_na(name) %>% 
  drop_na(fle_diet_prop) %>% 
  drop_na(cod_diet_prop) %>% 
  group_by(year, quarter, ices_rect) %>% 
  summarise(levin_cod = (1/(16-1)) * (((1/sum(cod_diet_prop^2))) - 1),
            levin_fle = (1/(16-1)) * (((1/sum(fle_diet_prop^2))) - 1)) %>% 
  ungroup()
#> Joining, by = c("name", "year", "quarter", "ices_rect")
#> left_join: added 3 columns (fle_diet, fle_diet_tot, fle_diet_prop)
#>            > rows only in x     432
#>            > rows only in y  (   32)
#>            > matched rows     1,264
#>            >                 =======
#>            > rows total       1,696
#> drop_na: no rows removed
#> drop_na: removed 448 rows (26%), 1,248 rows remaining
#> drop_na: no rows removed
#> group_by: 3 grouping variables (year, quarter, ices_rect)
#> summarise: now 78 rows and 5 columns, 2 group variables remaining (year, quarter)
#> ungroup: no grouping variables

# Merge the indicies  
ind <- left_join(levin, schoener)
#> Joining, by = c("year", "quarter", "ices_rect")
#> left_join: added one column (schoener)
#>            > rows only in x    0
#>            > rows only in y  ( 0)
#>            > matched rows     78
#>            >                 ====
#>            > rows total       78

# Summarise cod and flounder data by ices_rect then add to diet data
colnames(cod)
#>  [1] "FR_tot"                   "FR_sad"                  
#>  [3] "FR_spr"                   "saduria_entomon_tot"     
#>  [5] "tot_prey_biom"            "common_tot"              
#>  [7] "unique_pred_id"           "year"                    
#>  [9] "quarter"                  "time_period"             
#> [11] "quarter_fact"             "pred_weight_g"           
#> [13] "pred_cm"                  "predator_spec"           
#> [15] "predcod_density"          "predfle_density"         
#> [17] "predcod_density_sc"       "predfle_density_sc"      
#> [19] "depth"                    "X"                       
#> [21] "Y"                        "lat"                     
#> [23] "long"                     "ices_rect"               
#> [25] "cruise"                   "depth2_sc"               
#> [27] "saduria_entomon_per_mass" "tot_prey_biom_per_mass"  
#> [29] "depth_sc"
cod$year_rect_id <- paste(cod$year, cod$quarter, cod$ices_rect, sep = "_")

dens_sum <- cod %>% group_by(year_rect_id) %>%
  summarise(predfle_density_tot = sum(predfle_density),
            predcod_density_tot = sum(predcod_density)) %>% 
  ungroup() %>% 
  mutate(predfle_density_tot_sc = (predfle_density_tot - mean(predfle_density_tot)) / sd(predfle_density_tot),
         predcod_density_tot_sc = (predcod_density_tot - mean(predcod_density_tot)) / sd(predcod_density_tot))
#> group_by: one grouping variable (year_rect_id)
#> summarise: now 112 rows and 3 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'predfle_density_tot_sc' (double) with 110 unique values and 0% NA
#>         new variable 'predcod_density_tot_sc' (double) with 110 unique values and 0% NA

ind$year_rect_id <- paste(ind$year, ind$quarter, ind$ices_rect, sep = "_")

ind <- left_join(ind, dens_sum)
#> Joining, by = "year_rect_id"
#> left_join: added 4 columns (predfle_density_tot, predcod_density_tot, predfle_density_tot_sc, predcod_density_tot_sc)
#>            > rows only in x    0
#>            > rows only in y  (34)
#>            > matched rows     78
#>            >                 ====
#>            > rows total       78

# Summarise depth from the prediction grid then add to diet data
pred_grid2_sum <- pred_grid2 %>% 
  group_by(ices_rect) %>%
  summarise(mean_depth = mean(depth)) %>%
  mutate(depth_sc = (mean_depth - mean(mean_depth)) / sd(mean_depth)) %>% 
  ungroup()
#> group_by: one grouping variable (ices_rect)
#> summarise: now 61 rows and 2 columns, ungrouped
#> mutate: new variable 'depth_sc' (double) with 61 unique values and 0% NA
#> ungroup: no grouping variables

ind <- left_join(ind, dplyr::select(pred_grid2_sum, ices_rect, depth_sc, mean_depth))
#> Joining, by = "ices_rect"
#> left_join: added 2 columns (depth_sc, mean_depth)
#>            > rows only in x    0
#>            > rows only in y  (43)
#>            > matched rows     78
#>            >                 ====
#>            > rows total       78

# Plot the indicies
ggplot(ind) +
  geom_jitter(aes(factor(year), schoener),
              alpha = 0.8, width = 0.2, height = 0, size = 2)

ggplot(ind) +
  geom_jitter(aes(factor(year), levin_cod, color = "cod"),
              alpha = 0.8, width = 0.2, height = 0, size = 2) + 
  geom_jitter(aes(factor(year), levin_fle, color = "fle"),
              alpha = 0.8, width = 0.2, height = 0, size = 2) +
  scale_color_brewer(palette = "Dark2")


# Against response variables
# Schoener 
ggplot(ind, aes(x = predfle_density_tot + predcod_density_tot, y = schoener)) +
  geom_point() + stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter)

  
ggplot(ind, aes(x = predfle_density_tot, y = schoener)) + geom_point() + 
  stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter)


ggplot(ind, aes(x = predcod_density_tot, y = schoener)) + geom_point() +
  stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter)


# Now Levin
# Cod
ggplot(ind, aes(x = predfle_density_tot + predcod_density_tot, y = levin_cod)) +
  geom_point() + stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter) 

  
ggplot(ind, aes(x = predfle_density_tot, y = levin_cod)) + geom_point() +
  stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter) 

  
ggplot(ind, aes(x = predcod_density_tot, y = levin_cod)) + geom_point() + 
  stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter) 


# Flounder
ggplot(ind, aes(x = predfle_density_tot + predcod_density_tot, y = levin_fle)) +
  geom_point() + stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter) 

  
ggplot(ind, aes(x = predfle_density_tot, y = levin_fle)) + geom_point() +
  stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter) 

  
ggplot(ind, aes(x = predcod_density_tot, y = levin_fle)) + geom_point() + 
  stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter) 


# Add refs here or in text already! THis is complementary to prey weights in stomach
# Do cumsum plots
# Now split by size-classes, check Orior co author paper for which sizes to split
# Think about which model to fit to these indices?
# ---- brms with random rectangle and fixed year? Evaluate differences in intercept and plot covariates?

# ... Then when I have size-group? I get like 5 "species". Maybe do pca on those?
# Do that after resolving old comments, adding new diet data and correctly adding quarter and new cpue data!!!

# Quickly check data to determine which distribution to use
ind %>% 
  ungroup() %>% 
  count(schoener == 0) %>% 
  mutate(prop = n / sum(n))
#> ungroup: no grouping variables
#> count: now 2 rows and 2 columns, ungrouped
#> mutate: new variable 'prop' (double) with 2 unique values and 0% NA
#> # A tibble: 2 × 3
#>   `schoener == 0`     n   prop
#>   <lgl>           <int>  <dbl>
#> 1 FALSE              75 0.962 
#> 2 TRUE                3 0.0385

ind %>% 
  ungroup() %>% 
  count(levin_cod == 0) %>% 
  mutate(prop = n / sum(n))
#> ungroup: no grouping variables
#> count: now 2 rows and 2 columns, ungrouped
#> mutate: new variable 'prop' (double) with 2 unique values and 0% NA
#> # A tibble: 2 × 3
#>   `levin_cod == 0`     n   prop
#>   <lgl>            <int>  <dbl>
#> 1 FALSE               77 0.987 
#> 2 TRUE                 1 0.0128

ind %>% 
  ungroup() %>% 
  count(levin_fle == 0) %>% 
  mutate(prop = n / sum(n))
#> ungroup: no grouping variables
#> count: now 2 rows and 2 columns, ungrouped
#> mutate: new variable 'prop' (double) with 2 unique values and 0% NA
#> # A tibble: 2 × 3
#>   `levin_fle == 0`     n   prop
#>   <lgl>            <int>  <dbl>
#> 1 FALSE               75 0.962 
#> 2 TRUE                 3 0.0385

# Fit beta models, so few zeroes!
ind %>% arrange(schoener) %>% dplyr::select(schoener)
#> # A tibble: 78 × 1
#>    schoener
#>       <dbl>
#>  1 0       
#>  2 0       
#>  3 0       
#>  4 0.000490
#>  5 0.00127 
#>  6 0.00134 
#>  7 0.00149 
#>  8 0.00209 
#>  9 0.00240 
#> 10 0.00261 
#> # … with 68 more rows
ind %>% arrange(levin_cod) %>% dplyr::select(levin_cod)
#> # A tibble: 78 × 1
#>    levin_cod
#>        <dbl>
#>  1   0      
#>  2   0.00173
#>  3   0.00265
#>  4   0.00493
#>  5   0.00801
#>  6   0.0132 
#>  7   0.0163 
#>  8   0.0179 
#>  9   0.0191 
#> 10   0.0270 
#> # … with 68 more rows
ind %>% arrange(levin_fle) %>% dplyr::select(levin_fle)
#> # A tibble: 78 × 1
#>    levin_fle
#>        <dbl>
#>  1  0       
#>  2  0       
#>  3  0       
#>  4  0.000321
#>  5  0.000737
#>  6  0.00191 
#>  7  0.00206 
#>  8  0.00304 
#>  9  0.00779 
#> 10  0.00811 
#> # … with 68 more rows

ind <- ind %>% mutate(schoener2 = ifelse(schoener == 0, 0.0001, schoener),
                      levin_cod2 = ifelse(levin_cod == 0, 0.0001, levin_cod),
                      levin_fle2 = ifelse(levin_fle == 0, 0.0001, levin_fle),
                      year_f = as.factor(year),
                      quarter_f = as.factor(quarter),
                      ices_rect_f = as.factor(ices_rect))
#> mutate: new variable 'schoener2' (double) with 76 unique values and 0% NA
#>         new variable 'levin_cod2' (double) with 78 unique values and 0% NA
#>         new variable 'levin_fle2' (double) with 76 unique values and 0% NA
#>         new variable 'year_f' (factor) with 6 unique values and 0% NA
#>         new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#>         new variable 'ices_rect_f' (factor) with 18 unique values and 0% NA

Fit brms models with densities as covariates to diversity and overlap indices

# All covariates
m_schoen_beta_full <- brm(
  bf(schoener2 ~ 0 + year_f + quarter_f + depth_sc + predfle_density_tot_sc + predcod_density_tot_sc + (1|ices_rect_f),
     phi ~ 0 + year_f + quarter_f),
  data = ind, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
  #backend = "cmdstanr",
  chains = 4, iter = 4000, warmup = 1000, cores = 4, control = list(adapt_delta = 0.95))
#> Compiling Stan program...
#> Trying to compile a simple C file
#> Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
#> clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Users/maxlindmark/Library/R/4.0/library/Rcpp/include/"  -I"/Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/"  -I"/Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/unsupported"  -I"/Users/maxlindmark/Library/R/4.0/library/BH/include" -I"/Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/src/"  -I"/Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/"  -I"/Users/maxlindmark/Library/R/4.0/library/RcppParallel/include/"  -I"/Users/maxlindmark/Library/R/4.0/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
#> In file included from <built-in>:1:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Core:88:
#> /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
#> namespace Eigen {
#> ^
#> /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
#> namespace Eigen {
#>                ^
#>                ;
#> In file included from <built-in>:1:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
#> /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#> #include <complex>
#>          ^~~~~~~~~
#> 3 errors generated.
#> make: *** [foo.o] Error 1
#> Start sampling

plot(m_schoen_beta_full)

conditional_effects(m_schoen_beta_full)

loo_m_schoen_beta_full <- loo(m_schoen_beta_full, moment_match = TRUE)
#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
plot(loo_m_schoen_beta_full)


# Only density covariates
m_schoen_beta_dens <- brm(
  bf(schoener2 ~ 0 + year_f + predfle_density_tot_sc + predcod_density_tot_sc + (1|ices_rect_f),
     phi ~ 0 + year_f),
  data = ind, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
  #backend = "cmdstanr",
  chains = 4, iter = 4000, warmup = 1000, cores = 4, control = list(adapt_delta = 0.95))
#> Compiling Stan program...
#> Start sampling

plot(m_schoen_beta_dens)

conditional_effects(m_schoen_beta_dens)

loo_m_schoen_beta_dens <- loo(m_schoen_beta_dens, moment_match = TRUE)
#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
plot(loo_m_schoen_beta_dens)


# Simplest model
m_schoen_beta <- brm(
  bf(schoener2 ~ 0 + year_f + (1|ices_rect_f),
     phi ~ 0 + year_f),
  data = ind, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
  #backend = "cmdstanr",
  chains = 4, iter = 4000, warmup = 1000, cores = 4, control = list(adapt_delta = 0.95))
#> Compiling Stan program...
#> Start sampling

plot(m_schoen_beta)

conditional_effects(m_schoen_beta)

loo_m_schoen_beta <- loo(m_schoen_beta, moment_match = TRUE)
#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
plot(loo_m_schoen_beta)


loo_compare(loo_m_schoen_beta, loo_m_schoen_beta_dens, loo_m_schoen_beta_full)
#>                    elpd_diff se_diff
#> m_schoen_beta       0.0       0.0   
#> m_schoen_beta_dens -1.0       1.4   
#> m_schoen_beta_full -4.3       2.6

# Plot models (working with the m_schoen_beta_dens)
summary(m_schoen_beta_dens)
#>  Family: beta 
#>   Links: mu = logit; phi = log 
#> Formula: schoener2 ~ 0 + year_f + predfle_density_tot_sc + predcod_density_tot_sc + (1 | ices_rect_f) 
#>          phi ~ 0 + year_f
#>    Data: ind (Number of observations: 78) 
#>   Draws: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
#>          total post-warmup draws = 12000
#> 
#> Group-Level Effects: 
#> ~ices_rect_f (Number of levels: 18) 
#>               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept)     0.19      0.16     0.01     0.58 1.00     4846     5380
#> 
#> Population-Level Effects: 
#>                        Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> year_f2015                -2.74      0.46    -3.59    -1.75 1.00    10133
#> year_f2016                -2.91      0.30    -3.48    -2.29 1.00    11525
#> year_f2017                -2.57      0.36    -3.24    -1.83 1.00    11899
#> year_f2018                -2.14      0.51    -3.08    -1.08 1.00    10237
#> year_f2019                -1.58      0.51    -2.55    -0.55 1.00    12042
#> year_f2020                -1.46      0.42    -2.26    -0.60 1.00    11745
#> predfle_density_tot_sc    -0.25      0.16    -0.57     0.06 1.00    10524
#> predcod_density_tot_sc     0.18      0.15    -0.11     0.48 1.00    10331
#> phi_year_f2015             2.36      0.64     0.95     3.46 1.00     9842
#> phi_year_f2016             2.30      0.38     1.49     2.97 1.00    11640
#> phi_year_f2017             1.82      0.42     0.91     2.56 1.00    12310
#> phi_year_f2018             1.21      0.54     0.02     2.17 1.00    10486
#> phi_year_f2019             0.67      0.51    -0.40     1.58 1.00    11405
#> phi_year_f2020             1.46      0.50     0.40     2.35 1.00    12856
#>                        Tail_ESS
#> year_f2015                 7731
#> year_f2016                 7531
#> year_f2017                 8457
#> year_f2018                 7513
#> year_f2019                 8082
#> year_f2020                 8680
#> predfle_density_tot_sc     9258
#> predcod_density_tot_sc     8973
#> phi_year_f2015             7032
#> phi_year_f2016             7508
#> phi_year_f2017             7677
#> phi_year_f2018             7340
#> phi_year_f2019             8729
#> phi_year_f2020             8334
#> 
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

# Working with the posterior
posterior_beta <- m_schoen_beta_dens %>% 
  gather_draws(`b_.*`, regex = TRUE) %>% 
  mutate(component = ifelse(str_detect(.variable, "phi_"), "Precision", "Mean"),
         intercept = str_detect(.variable, "Intercept"))

ggplot(posterior_beta, aes(x = .value, y = fct_rev(.variable), fill = component)) +
  geom_vline(xintercept = 0) +
  stat_halfeye(aes(slab_alpha = intercept), alpha = 0.5,
               .width = c(0.8, 0.95), point_interval = "median_hdi") +
  scale_fill_brewer(palette = "Dark2") +
  scale_slab_alpha_discrete(range = c(1, 0.4)) +
  guides(fill = "none", slab_alpha = "none") +
  labs(x = "Coefficient", y = "Variable") +
  facet_wrap(vars(component), ncol = 1, scales = "free_y") +
  theme_plot() + 
  NULL


# Selected parameters
# https://www.sciencedirect.com/topics/mathematics/logit-scale
posterior_beta %>%
  filter(component == "Mean" & .variable %in% c("b_predcod_density_tot_sc", "b_predfle_density_tot_sc")) %>% 
  ggplot(., aes(x = plogis(.value), y = .variable, fill = .variable)) +
  geom_vline(xintercept = 0) +
  stat_halfeye(alpha = 0.5, .width = c(0.8, 0.95)) +
  scale_fill_brewer(palette = "Dark2") +
  labs(x = "Coefficient", y = "Variable") +
  theme_plot() + 
  guides(fill = FALSE) + 
  NULL
#> filter (grouped): removed 144,000 rows (86%), 24,000 rows remaining
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.


# Prediction! (non dens model)
beta_bayes_pred_yr <- m_schoen_beta %>% 
  epred_draws(newdata = tibble(year_f = as.factor(c(2015:2020))),
              re_formula = NA)

ggplot(beta_bayes_pred_yr, aes(x = .epred, y = year_f, color = year_f, fill = year_f)) +
  stat_halfeye(.width = c(0.8, 0.95), alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(fill = FALSE, color = FALSE) +
  coord_flip(xlim = c(0, 0.4)) + 
  labs(x = "Predicted Schoener's overlap index", y = NULL) +
  theme_plot()
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.


beta_bayes_pred_re <- m_schoen_beta %>% 
  epred_draws(newdata = tibble(expand.grid(year_f = as.factor(c(2015:2020)),
                                            ices_rect_f = as.factor(unique(ind$ices_rect)))))

ggplot(beta_bayes_pred_re, aes(x = .epred, y = ices_rect_f, color = ices_rect_f, fill = ices_rect_f)) +
  stat_halfeye(.width = c(0.8, 0.95), alpha = 0.7) +
  scale_fill_viridis(discrete = TRUE) +
  scale_color_viridis(discrete = TRUE) +
  guides(fill = FALSE, color = FALSE) +
  coord_cartesian(xlim = c(0, 0.4)) + 
  #coord_flip(xlim = c(0, 0.4)) + 
  labs(x = "Predicted Schoener's overlap index", y = NULL) +
  theme_plot()
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.

# Cod Levin index: all covariates
m_levin_beta_full_cod <- brm(
  bf(levin_cod2 ~ 0 + year_f + quarter_f + depth_sc + predfle_density_tot_sc + predcod_density_tot_sc + (1|ices_rect_f),
     phi ~ 0 + year_f + quarter_f),
  data = ind, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
  #backend = "cmdstanr",
  chains = 2, iter = 2000, warmup = 1000, cores = 2, control = list(adapt_delta = 0.99))
#> Compiling Stan program...
#> Start sampling

plot(m_levin_beta_full_cod)

conditional_effects(m_levin_beta_full_cod)

loo_m_schoen_beta_full <- loo(m_levin_beta_full_cod, moment_match = TRUE)
#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
plot(m_levin_beta_full_cod)


m_levin_beta_full_fle <- brm(
  bf(levin_fle2 ~ 0 + year_f + quarter_f + depth_sc + predfle_density_tot_sc + predcod_density_tot_sc + (1|ices_rect_f),
     phi ~ 0 + year_f + quarter_f),
  data = ind, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
  #backend = "cmdstanr",
  chains = 2, iter = 2000, warmup = 1000, cores = 2, control = list(adapt_delta = 0.99))
#> Compiling Stan program...
#> Start sampling

plot(m_levin_beta_full_fle)

conditional_effects(m_levin_beta_full_fle)

loo_m_schoen_beta_full <- loo(m_levin_beta_full_fle, moment_match = TRUE)
#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
plot(m_levin_beta_full_fle)

Fit sdmTMB models with densities as covariates to stomach content data. First make mesh

# Cod 
pcod_spde <- make_mesh(cod, c("X", "Y"), n_knots = 150, type = "kmeans", seed = 42)
plot(pcod_spde)


pfle_spde <- make_mesh(fle, c("X", "Y"), n_knots = 70, type = "kmeans", seed = 42)
plot(pfle_spde)

Fit models of total content

# Model total prey biomass 
# Cod
mcod <- sdmTMB(
  data = cod, 
  formula = tot_prey_biom_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + predcod_density_sc + s(depth_sc),
  time = "year", fields = "IID", spatial_only = FALSE, spde = pcod_spde, family = tweedie())
#> Registered S3 methods overwritten by 'lme4':
#>   method                          from
#>   cooks.distance.influence.merMod car 
#>   influence.merMod                car 
#>   dfbeta.influence.merMod         car 
#>   dfbetas.influence.merMod        car
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k` since the degree of wiggliness
#> is determined by the data.

tidy(mcod)
#>                  term    estimate  std.error
#> 1            quarter1 -4.81710803 0.18724912
#> 2            quarter4 -4.66825835 0.16127572
#> 3 as.factor(year)2016  0.24009064 0.19659901
#> 4 as.factor(year)2017  0.36443934 0.19436610
#> 5 as.factor(year)2018  0.42481751 0.21019866
#> 6 as.factor(year)2019  0.23688896 0.23211806
#> 7 as.factor(year)2020  0.47897776 0.23197072
#> 8  predfle_density_sc -0.02663818 0.02644890
#> 9  predcod_density_sc  0.09099320 0.03395342
summary(mcod)
#> Spatial model fit by ML ['sdmTMB']
#> Formula: tot_prey_biom_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + 
#> Formula:     predcod_density_sc + s(depth_sc)
#> Time column: "year"
#> SPDE: pcod_spde
#> Data: cod
#> Family: tweedie(link = 'log')
#>                     coef.est coef.se
#> quarter1               -4.82    0.19
#> quarter4               -4.67    0.16
#> as.factor(year)2016     0.24    0.20
#> as.factor(year)2017     0.36    0.19
#> as.factor(year)2018     0.42    0.21
#> as.factor(year)2019     0.24    0.23
#> as.factor(year)2020     0.48    0.23
#> predfle_density_sc     -0.03    0.03
#> predcod_density_sc      0.09    0.03
#> 
#> Dispersion parameter: 0.52
#> Tweedie p: 1.74
#> Matern range: 9.74
#> Spatial SD: 0.23
#> Spatiotemporal SD: 0.68
#> ML criterion at convergence: -9961.622
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.

cod$resids_mcod <- residuals(mcod) # randomized quantile residuals
qqnorm(cod$resids_mcod); abline(a = 0, b = 1)


# Flounder
mfle <- sdmTMB(
  data = fle, 
  formula = tot_prey_biom_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + predcod_density_sc + s(depth_sc),
  time = "year", fields = "IID", spatial_only = FALSE, spde = pfle_spde, family = tweedie())
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k` since the degree of wiggliness
#> is determined by the data.

tidy(mfle)
#>                  term     estimate  std.error
#> 1            quarter1 -5.038155329 0.22043580
#> 2            quarter4 -4.544660978 0.18833024
#> 3 as.factor(year)2016 -0.020008080 0.24180860
#> 4 as.factor(year)2017  0.378980202 0.23768869
#> 5 as.factor(year)2018  0.593044431 0.28817619
#> 6 as.factor(year)2019  0.584962410 0.25846307
#> 7 as.factor(year)2020  0.995539716 0.25328069
#> 8  predfle_density_sc -0.043565572 0.03287689
#> 9  predcod_density_sc  0.002175145 0.03343344
summary(mfle)
#> Spatial model fit by ML ['sdmTMB']
#> Formula: tot_prey_biom_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + 
#> Formula:     predcod_density_sc + s(depth_sc)
#> Time column: "year"
#> SPDE: pfle_spde
#> Data: fle
#> Family: tweedie(link = 'log')
#>                     coef.est coef.se
#> quarter1               -5.04    0.22
#> quarter4               -4.54    0.19
#> as.factor(year)2016    -0.02    0.24
#> as.factor(year)2017     0.38    0.24
#> as.factor(year)2018     0.59    0.29
#> as.factor(year)2019     0.58    0.26
#> as.factor(year)2020     1.00    0.25
#> predfle_density_sc     -0.04    0.03
#> predcod_density_sc      0.00    0.03
#> 
#> Dispersion parameter: 0.12
#> Tweedie p: 1.50
#> Matern range: 29.17
#> Spatial SD: 0.30
#> Spatiotemporal SD: 0.54
#> ML criterion at convergence: -4721.832
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.

fle$resids_mfle <- residuals(mfle) # randomized quantile residuals
qqnorm(fle$resids_mfle); abline(a = 0, b = 1)

Fit models of saduria content

# Model Saduria contents 
# Cod
mcodsad <- sdmTMB(
  data = cod, 
  formula = saduria_entomon_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + predcod_density_sc + s(depth_sc),
  time = "year", fields = "IID", spatial_only = FALSE, spde = pcod_spde, family = tweedie())
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k` since the degree of wiggliness
#> is determined by the data.

tidy(mcodsad)
#>                  term    estimate std.error
#> 1            quarter1 -9.02384114 0.7872885
#> 2            quarter4 -9.02682697 0.7395658
#> 3 as.factor(year)2016 -1.07056118 0.7875507
#> 4 as.factor(year)2017 -1.45889570 0.7842763
#> 5 as.factor(year)2018 -0.21692129 0.8426281
#> 6 as.factor(year)2019 -1.04243883 0.9148954
#> 7 as.factor(year)2020  0.06204430 0.8954815
#> 8  predfle_density_sc -0.10384084 0.1238834
#> 9  predcod_density_sc -0.04374943 0.1062921
summary(mcodsad)
#> Spatial model fit by ML ['sdmTMB']
#> Formula: saduria_entomon_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + 
#> Formula:     predcod_density_sc + s(depth_sc)
#> Time column: "year"
#> SPDE: pcod_spde
#> Data: cod
#> Family: tweedie(link = 'log')
#>                     coef.est coef.se
#> quarter1               -9.02    0.79
#> quarter4               -9.03    0.74
#> as.factor(year)2016    -1.07    0.79
#> as.factor(year)2017    -1.46    0.78
#> as.factor(year)2018    -0.22    0.84
#> as.factor(year)2019    -1.04    0.91
#> as.factor(year)2020     0.06    0.90
#> predfle_density_sc     -0.10    0.12
#> predcod_density_sc     -0.04    0.11
#> 
#> Dispersion parameter: 0.17
#> Tweedie p: 1.49
#> Matern range: 41.79
#> Spatial SD: 2.00
#> Spatiotemporal SD: 1.53
#> ML criterion at convergence: -1022.716
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.

cod$resids_mcodsad <- residuals(mcodsad) # randomized quantile residuals
qqnorm(cod$resids_mcodsad); abline(a = 0, b = 1)


# Flounder
mflesad <- sdmTMB(
  data = fle, 
  formula = saduria_entomon_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + predcod_density_sc + s(depth_sc),
  time = "year", fields = "IID", spatial_only = FALSE, spde = pfle_spde, family = tweedie())
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k` since the degree of wiggliness
#> is determined by the data.

tidy(mflesad)
#>                  term    estimate  std.error
#> 1            quarter1 -8.99399061 0.98263924
#> 2            quarter4 -7.87852100 0.94515317
#> 3 as.factor(year)2016  1.00546603 0.80258800
#> 4 as.factor(year)2017  0.99471876 0.79567385
#> 5 as.factor(year)2018  1.89586027 0.91527579
#> 6 as.factor(year)2019  0.52603927 0.85483508
#> 7 as.factor(year)2020  2.15101984 0.82404480
#> 8  predfle_density_sc  0.02410729 0.07481394
#> 9  predcod_density_sc -0.17705506 0.07777392
summary(mflesad)
#> Spatial model fit by ML ['sdmTMB']
#> Formula: saduria_entomon_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + 
#> Formula:     predcod_density_sc + s(depth_sc)
#> Time column: "year"
#> SPDE: pfle_spde
#> Data: fle
#> Family: tweedie(link = 'log')
#>                     coef.est coef.se
#> quarter1               -8.99    0.98
#> quarter4               -7.88    0.95
#> as.factor(year)2016     1.01    0.80
#> as.factor(year)2017     0.99    0.80
#> as.factor(year)2018     1.90    0.92
#> as.factor(year)2019     0.53    0.85
#> as.factor(year)2020     2.15    0.82
#> predfle_density_sc      0.02    0.07
#> predcod_density_sc     -0.18    0.08
#> 
#> Dispersion parameter: 0.14
#> Tweedie p: 1.47
#> Matern range: 77.53
#> Spatial SD: 2.11
#> Spatiotemporal SD: 1.18
#> ML criterion at convergence: -1467.275
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.

fle$resids_mflesad <- residuals(mflesad) # randomized quantile residuals
qqnorm(fle$resids_mflesad); abline(a = 0, b = 1)

Calculate and plot marginal effects

Cod

# Flounder density
nd_cod_fle <- data.frame(predfle_density_sc = seq(min(cod$predfle_density_sc),
                                                  max(cod$predfle_density_sc), length.out = 100))

nd_cod_fle$year <- 2018L
nd_cod_fle$predcod_density_sc <- 0
nd_cod_fle$depth_sc <- 0
nd_cod_fle$quarter <- factor(4)

sad_pred_cod_fle <- predict(mcodsad, newdata = nd_cod_fle, se_fit = TRUE, re_form = NA)
tot_pred_cod_fle <- predict(mcod, newdata = nd_cod_fle, se_fit = TRUE, re_form = NA)

# Cod density
nd_cod_cod <- data.frame(predcod_density_sc = seq(min(cod$predcod_density_sc),
                                                  max(cod$predcod_density_sc), length.out = 100))

nd_cod_cod$year <- 2018L
nd_cod_cod$predfle_density_sc <- 0
nd_cod_cod$depth_sc <- 0
nd_cod_cod$quarter <- factor(4)

sad_pred_cod_cod <- predict(mcodsad, newdata = nd_cod_cod, se_fit = TRUE, re_form = NA)
tot_pred_cod_cod <- predict(mcod, newdata = nd_cod_cod, se_fit = TRUE, re_form = NA)

Flounder

# Flounder density
nd_fle_fle <- data.frame(predfle_density_sc = seq(min(fle$predfle_density_sc),
                                                  max(fle$predfle_density_sc), length.out = 100))

nd_fle_fle$year <- 2018L
nd_fle_fle$predcod_density_sc <- 0
nd_fle_fle$depth_sc <- 0
nd_fle_fle$quarter <- factor(4)

sad_pred_fle_fle <- predict(mflesad, newdata = nd_fle_fle, se_fit = TRUE, re_form = NA)
tot_pred_fle_fle <- predict(mfle, newdata = nd_fle_fle, se_fit = TRUE, re_form = NA)

# Cod density
nd_fle_cod <- data.frame(predcod_density_sc = seq(min(fle$predcod_density_sc),
                                                  max(fle$predcod_density_sc), length.out = 100))

nd_fle_cod$year <- 2018L
nd_fle_cod$predfle_density_sc <- 0
nd_fle_cod$depth_sc <- 0
nd_fle_cod$quarter <- factor(4)

sad_pred_fle_cod <- predict(mflesad, newdata = nd_fle_cod, se_fit = TRUE, re_form = NA)
tot_pred_fle_cod <- predict(mfle, newdata = nd_fle_cod, se_fit = TRUE, re_form = NA)

Plot marginal effects

# Saduria models
p1 <- ggplot(sad_pred_fle_fle, aes(predfle_density_sc, exp(est),
  ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4) +
  #coord_cartesian(ylim = c(0, 0.015), expand = 0) + 
  labs(y = "Saduria in flounder stomach [g/g]", x = "")

p2 <- ggplot(sad_pred_fle_cod, aes(predcod_density_sc, exp(est),
  ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4) + 
  coord_cartesian(ylim = c(0, 0.015), expand = 0) + 
  labs(y = "", x = "")

p3 <- ggplot(sad_pred_cod_fle, aes(predfle_density_sc, exp(est),
  ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4) +
  labs(y = "Saduria in cod stomach [g/g]", x = "Scaled flounder density")

p4 <- ggplot(sad_pred_cod_cod, aes(predcod_density_sc, exp(est),
  ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4) + 
  labs(y = "", x = "Scaled cod density")

p1 + p2 + p3 + p4 + plot_layout(ncol = 2)


ggsave("figures/marginal_effects_saduria.png", width = 6.5, height = 6.5, dpi = 600)

# Total stomach content models
p5 <- ggplot(tot_pred_fle_fle, aes(predfle_density_sc, exp(est),
  ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4) +
  coord_cartesian(ylim = c(0, 0.033), expand = 0) + 
  labs(y = "Prey biomass in flounder stomach [g/g]", x = "")

p6 <- ggplot(tot_pred_fle_cod, aes(predcod_density_sc, exp(est),
  ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4) + 
  coord_cartesian(ylim = c(0, 0.033), expand = 0) + 
  labs(y = "", x = "")

p7 <- ggplot(tot_pred_cod_fle, aes(predfle_density_sc, exp(est),
  ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4) +
  coord_cartesian(ylim = c(0, 0.028), expand = 0) + 
  labs(y = "Prey biomass in cod stomach [g/g]", x = "Scaled flounder density")

p8 <- ggplot(tot_pred_cod_cod, aes(predcod_density_sc, exp(est),
  ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
  geom_line() + geom_ribbon(alpha = 0.4) + 
  coord_cartesian(ylim = c(0, 0.028), expand = 0) + 
  labs(y = "", x = "Scaled cod density")

p5 + p6 + p7 + p8 + plot_layout(ncol = 2)


ggsave("figures/marginal_effects_total.png", width = 6.5, height = 6.5, dpi = 600)
knitr::knit_exit()